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ABSTRACT 

We develop a framework for constructing mixed multiscale finite volume methods for 
elliptic equations with multiple scales arising from flows in porous media. Some of the 
methods developed using the framework are already known [2D] ; others are new. New insight 
is gained for the known methods and extra flexibility is provided by the new methods. We 
give as an example a mixed MsFV on uniform mesh in 2-D. This method uses novel multiscale 
velocity basis functions that are suited for using global information, which is often needed 
to improve the accuracy of the multiscale simulations in the case of continuum scales with 
strong non-local features. The method efficiently captures the small effects on a coarse 
grid. We analyze the new mixed MsFV and apply it to solve two-phase flow equations in 
heterogeneous porous media. Numerical examples demonstrate the accuracy and efficiency 
of the proposed method for modeling the flows in porous media with non-separable and 
separable scales. 

1 Introduction 

Subsurface flows are often affected by heterogeneities in a wide range of length scales. This 
causes significant challenges for subsurface flow modeling. Geological characterizations that 
capture these effects are typically developed at scales that are too fine for direct flow simu- 
lations. Usually, upscaled or multiscale models are employed for such systems. In upscaling 
methods, the original model is coarsened by numerically homogenizing parameters (e.g., 
permeability). The simulation is performed using the coarsened model, which may differ 
from the underlying fine-scale model. In multiscale methods, the fine-scale information is 
carried throughout the simulation and the coarse-scale equations are generally not expressed 
analytically, but rather formed and solved numerically. 

Various numerical multiscale approaches for flows in porous media have been developed 
during the past decade. A multiscale finite element method (MsFEM) was introduced in 
[IS] and takes its origin from the pioneering work [5]. Its main idea is to incorporate the 
small-scale information into finite element basis functions and capture their effect on the 
large scales via finite element computations. The MsFEM in [TS] shares some similarities 
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with a number of multiscale numerical methods, such as residual free bubbles [6], variational 
multiscale method [19], two-scale conservative subgrid approaches [3], heterogeneous multi- 
scale method [15] and multiscale discontinuous Galerkin method [31] . Chen and Hou have 
applied the MsFEM idea in combination with a mixed finite element formulation to propose 
a mixed MsFEM [8]. Recently, Arbogast et al. [I] used domain decomposition approach 
and variational mixed formulation to develop a multiscale mortar mixed MsFEM. Jenny et 
al. [20] have used the ideas in [18] and finite volume framework to design a multiscale fi- 
nite volume method (MsFV). The MsFV and its variants have proved successful in reservoir 
simulations. 

Here we develop a framework for constructing mixed MsFV methods, which uses ideas 
from the mixed finite volume methods [2U [251 I2E] , multi-point flux approximations (MPFA) 
[21 [16], and mixed MsFEM. The mixed MsFV are mass conservative methods, which is 
an important property of the discretizations used in subsurface flow simulations (see [TT] 
for related discussion). The important feature of the mixed finite volume methods is the 
direct approximation of the velocity, that is, specially constructed discrete spaces are used 
to approximate the velocity unknowns. We propose a novel way to construct multiscale 
velocity basis functions that are well suited for parallel computation. Mixed MsFEM reduces 
the system of coupled equations for pressure and velocity to a system only for the pressure. 
However, the reduction process is computational expensive and has some restrictions when 
the global mass matrix in mixed MsFEM is large. In the mixed MsFV, we compute the 
inverse of each local mass matrix instead of global mass matrix and get effective coarse- 
scale transmissibilities. This computation is cheap and efficient. In the MsFV proposed 
in [20] , two sets of multiscale basis functions are computed: the first set of basis functions 
is to approximate pressure and the second set of basis functions is required to construct a 
conservative fine-scale velocity field. Only one set of multiscale basis functions is constructed 
in the mixed MsFV and the span of the basis functions are to approximate the velocity. 
Piecewise constant is used for pressure basis in the mixed MsFV. Hence the computation 
for basis functions in the mixed MsFV is less expensive than the MsFV. To the best of our 
knowledge, the mixed MsFV is a new numerical multiscale method. 

Boundary effect is a great issue in many multiscale methods (e.g., [T81 IE])- When we 
construct the multiscale basis functions in the mixed MsFV, we can use constant bound- 
ary condition for basis equations and obtain local multiscale basis. We find that the local 
multiscale basis in the mixed MsFV has the similar merit to the multiscale basis using over- 
sampling technique developed in [18] and they are able to reduce the boundary effect greatly. 
The mixed MsFV using local multiscale basis works well for most multiscale problems and 
particularly for the case of separable scales (e.g., periodic media). Furthermore, we can 
employ global information for the multiscale basis functions of the mixed MsFV. The global 
information usually represents long range features of flows and is used to construct multiscale 
basis functions. The global information is needed in the case of strong non-separable scales 
and renders much better accuracy than local multiscale approaches pp. 

The proposed mixed MsFV in some extend inherits the advantages of the mixed MsFEM 
and MsFV and alleviates the drawbacks of them without increase of the computational cost. 
Moreover, this method and its generalizations are well suited for computation on unstruc- 
tured grids and can be easily incorporated in production reservoir simulators. For example, 
to construct the velocity basis no geometric information is necessary. The pressure basis can 
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be computed using only the local matrices if it consists of discrete harmonic functions as 
demonstrated in [TJ]. Detailed description is provided in [27]. 

The rest of the paper is organized as follows. Section 2 is devoted to formulating a 
standard mixed finite volume method for a model elliptic equation. In Section 3, we apply 
the methodology of the mixed finite volume method to develop a mixed MsFV method. 
Here we design a new multiscale velocity basis function, analyze the proposed mixed MsFV, 
and address some computational issues. In Section 4, we generalize the approach applied in 
Section 3 to derive our first mixed MsFV method and show how given a mixed FV method 
we can derived a corresponding mixed MsFV. Several examples are given to demonstrate 
how the framework can be used to develop new mixed MsFV methods. In Section 5, we 
apply a mixed MsFV to incompressible two-phase flows in porous media with continuum 
scales and separable scales. Finally, some comments and conclusions are made. 



2 Mixed finite volume method formulation for a model 
problem 

We first define notations for function spaces used in the paper. 

L\Q) = {f(x)\ [ \f(x)\ 2 dx<oo}, 
Jn 

H\n) = {f(x)\f(x)eL 2 (n) and V/(x) G [L 2 (f2)] d }, 
H(div,tt) = {f(x)\f(x) G {L 2 (tt)] d and div(/(x)) G L 2 (f2)}. 

We consider the following model elliptic equation, 

— div(/c(x)Vp) = f(x) in Q, 
k{x)Vp ■ n = on dfl, 

pdx = 0, 



(2.1; 



where Q is a domain in M. d , d = 2 or 3 and / G L 2 (Q). Equation (12. ip is used to model many 
physical processes, for example, fluid flow in porous media. The coefficient k(x) represents 
the permeability and is often heterogeneous. Here p represents pressure. 

We define velocity u(x) = —k(x)Vp. To simplify presentation, we shall not write the 
spatial variables x for functions when no ambiguity occurs. Then (12. ip can be rewritten as 
a first order system 



k~ l u + Vp = 
div(n) = /. 



(2.2) 



The weak mixed formulation for (12.21) is: 
Find {u,p} G U x V such that 



Jn 

/ div(u)qdx = / fqdx Vg G <2, 
Jn Jn 



(2.3) 
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where U = H{div, Q), V = H\Q), V = (L 2 (Q)) d , and Q = L 2 (Q). 

There are different ways to construct a discretization of equation ( 12. 3p . One can try 
to find the approximations such that {uh,Ph} £ H(div,Q) x [30]. Unfortunately 

these discretizations are computationally expensive and are only applicable on very restric- 
tive grids. Frequently dual mixed finite element methods are used [TJ [29] to construct the 
discretization {uh,Ph} £ H(div,Q) x L 2 (Q). We will consider a class of methods that are 
related to the primal mixed finite element methods [29], i.e., the discretizations we seek are 
{uh,Ph} £ Uh x Vh, Uh C (L 2 (f2)) , Ph. C H l {VL) with additional conservation enforced on 
particular volumes. These methods are close related to the standard cell-centered finite vol- 
ume method and several multi-point flux approximation (MPFA) methods that generalize 
it. We will refer such approximations as mixed finite volume methods. 

We assume that two grids are defined: primary grid Th and dual grid Vh- Usually the 
primary grid is used to approximate the scaler variable p, and the dual mesh is used to 
construct the discretization of the velocity. Different examples of primary and dual grids 
can be found in j2j [TBI [2H EHl EE]. Consider the discrete problem: 

Find {uh,Ph} x Vh such that 



A particular mixed finite volume method will be fully described if we define the grids, the 
approximation spaces, and the operators V^, and div^. The operator div^ can be defined in 
the following way: 



with Qh the space of piece-wise constants on the volumes V from the primary grid and E is 
an edge in 2-D or face in 3-D in dV. Here rig denotes the outward normal vector to E. 

In the paper, we focus on 2-D case. The 3-D case is a straightforward extension of 2-D 
case. 

We give an example of mixed finite volume methods as following. 

Example 1 (Mixed FV 1). We assume that the primary grid Th consists of rectangles, and 
the dual grid T>h is also rectangular, but with vertexes the cell centers. The discrete space for 
the scalar variable, Vh = Qh, consists of piecewise constants and is defined on the primary 
grid, The approximation space for the vector variable, Uh = Vh, is the space of piece-wise 
constants vector constants with continuous normal components and is defined on the dual 
mesh. Consider one dual cell D = D { U Dj U D k U D t (see Figure \2.1\) . The four functions, 
£ik, £ji and eui are defined with the relations 



where pq and rt can be any element of the set Id = {ij, ik, jl, jk}. It is easy to see that 
{ e pg}y PQ £ Id ar G linearly independent and therefore form a basis ofUh\D- The degrees of 
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Figure 2.1: A dual cell D 



freedom are the integrals of the flux, i. e, the numbers v pq = J, v • n pq dx pq G Id- Then for 
any v e U h , v\ D = J2 pqe i D 

v pq e pq . The operator is given by: 



/ V h p h ■ v h dx = ^ / Vh ' Ue ^ e dx ' 



(2.6) 



where 



[p] E = lim(p(x + tn E ) - p(x - tn E )) 



t-fO 
t>0 



and the direction of n E is from left to right or from bottom to top in 2-D. Note that from the 
first equation of (12. 4p we can express 



Mu = p, 



(2.7) 



with u T = [Uij, Uik, Uj h u kl ] T , p = \pi~Pj,Pi-pk,Pj -pi,Pk-pi] T , and M a 4 x 4 matrix. We 
solve for u and plug the result in the second equation of (12. 4p to get the final discretization. 



For Example 1, We can rewrite equation (12. 5 p as 

/ Uh ' n El a h}E = ~ fqhdx. 



VeT h EeV ' 

which shows that the matrix of the discretization is symmetric. We note that Vh in Example 
1 is not subspace of H l (Q). Our approximation of Vp is nonconforming. 

The Example 1 is a straightforward generalization of the standard cell-centered finite 
volume method on structured rectangular meshes [25] . In fact, if the coefficient k(x) in (12.11) 
is a scalar function, then this method coincides with the standard cell-centered finite volume 
method. 

We need the matrix M in equation (I2.7P to have the appropriate dimensions and to 
be invertible in order to have a well defined discretization in Example 1. We state these 
condition for future reference: 

1. dim(Vh) + dim{Uh) = dim(Qh) + dim(Vh)] 

2. Matrix M is invertible. 

We will follow Example 1 to derive a new mixed MsFV method in the next section. 
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3 A new mixed multiscale finite volume method 



In order to describe the multiscale finite element method, we assume that the grids Th and 
T>h are coarse grids and that there exist an underlying fine grid containing the fine-scale 
information. Figure 13.21 depicts the rectangle primary coarse grid and dual coarse grid. The 
velocity u is discretized on the interfaces Eh of primary grid, e.g., Eij in Figure |3T2| and the 
pressure p is discretized on cell-centers of the primary mesh / the vertices of the dual grid, 
e.g, vertex k, I in Figure I3T21 

Define the discrete space for pressure to be Vh = Qh, the space of piecewise constants 
on Th, and for velocity to be Uh — V#, a multiscale finite element space with continuous 
normal on Eh that will be defined in Subsection 13. 1[ Then the mixed multiscale finite volume 
formulation for Equation (12. 3p reads: 

Find {uh,Ph} G Uh X Vh such that 



/ k 1 u H ■ v H dx + / v h ■ n E [pH\Edx = \/v H G U H , 

-* n EaEu e 



E 



(3i 



u H ■ nq H dx = / fq H dx Wq H G Vh, 
dKi Jn 



where \ph]e is the jump of pn across the interface E and defined in Example 1. It is clear 
from equation (13. 8p that the operators div# and Vh are defines as follows: 

(3.9) 



and 



/ VhPh -v H dx= ^2 / v h ■ n E [ph]e dx, 
/ div h (u h )q h dx = ^ ^ / u h -n E qhds (3.10) 



VeT h EedV 

2fr>\ n „J r- r2 



Clearly uh G L (Q) and pn G L (Q). By a similar argument as in Example 1, the discrete 
system of (13. 8p is symmetric. 



3.1 A new multiscale velocity basis function 



In order to complete the derivation of the method from (13.81) . we need to define velocity 
basis functions. In this subsection we design multiscale velocity basis functions associated 
with interfaces of Eh- 

Let e Eh be any interface and Qij (green part in Figure 13.21) be an open set bounded 



by edges e 



-ij, 



and 



interface as following: 



We construct a multiscale basis function <f>ij associated to the 



' — d\w{kV 4>ij) = in fi^, 

v(x)-n 



-kV<pi 



n 



t — i \ ^ on 

J e l v{x)-ndx ij 

ip 

v(x)-n 



J e l v(x)-ndx 





on e 



on e\jUe\j, 



(3.11) 
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Figure 3.2: Rectangle primary grid and dual grid 

where v(x) is a vector function and has some options depending on the multiscale features 
(e.g., separable scales or non-separable scales). We will address the options for v(x) later. 
Here n is the unit normal vector pointing out of Qij. We would like to note that the basis 
equation ( 13. lip is defined for a horizontal flux. By switching the no-flow boundary condition 
and flow boundary condition, we can similarly define the basis equation presenting a vertical 
flux. We define the velocity basis function ipij = —kV4>ij and the finite dimension space for 
velocity as 

Figure 13.31 depicts the vector fields of velocity basis functions (horizontal flux) for homoge- 
neous permeability and heterogeneous permeability (SPE 10, layer 85), respectively. The 
figure confirms that the multiscale basis defined in (13.111) reflects the properties of the me- 
dia/permeability. The multiscale basis functions are pre-computed and suited for parallel 
computation. 

For many multiscale problems, in particular, the problems with separable scales, we can 
simply take v(x) = (1, 1). No global information is used in this case, we call this case the 
local multiscale method. If the permeability k has strong long range features (e.g., highly 
channelized), we can use some global information v (x) related to the features of u, e.g., the 
single-phase velocity at time zero, to construct the velocity multiscale basis functions p]. 
In general, global information used to construct the basis functions in highly heterogeneous 
permeability often yields much better approximation than a constant boundary condition. 
Many numerical studies [U show the global information is very helpful to improve accu- 
racy when the permeability k has distinguished long range features. We note that one can 
use time-dependent global fields to construct multiscale basis functions. This is helpful for 
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Figure 3.3: Vector field of a velocity basis function for homogeneous permeability (left) and 
vector field of a velocity basis function for heterogeneous permeability (right) 



compressible flow simulations. 

If multiple global fields v n (x) {n 



N) are used to build basis functions, then each 



interface corresponds to multiple multiscale basis functions 



In 



-div(W0") = in 



■WC • n = { 



J i v n (x)-ndx 

v n (x)-n 
f i v n (x)-ndx 





on e 



on e- 



N), which solve 



(3.12) 



on eljUelj. 



We note that the multiple global fields v n (x) can be associated to some representative real- 
izations in the setting of stochastic two-phase flows [22]. 

Remark 3.1. By the result of Owhadi and Zhang [28], we can use d (d = dim(Q)) global 
fields in (13.121) . From the result, the global field v n (x) = — kVpn {n — 1, • ■ • , d) are the 
solutions of the elliptic equations 



-div(kVp n ) = in Q 

p n = x n on dfl, 



(3.13) 



where x — (x±, ■ ■ ■ , Xd). 

Following the ideas in [21], the global information can be computed on an intermediate 
coarse grid using upscaling techniques. This will reduce the computation for the global 
information. 

Because the source term of the basis equation (13. lip is zero, the flux conservation implies 
the following proposition. 

Proposition 3.1. Let e £h be any interface and ipij be the corresponding velocity basis 
function. Then 

ipij ■ n Eii dx = 1. 



Moreover, allipij are linear independent, i.e., they form a set of finite element basis functions. 



S 



Remark 3.2. If k is a constant in Qij and no global information is used, then ipij and ipn~ are 
orthogonal each other, i.e., the velocity basis function producing horizontal flux is orthogonal 
to the velocity basis function producing vertical flux. In fact, we can show that (see [25] ) 



, hi = { pferd.0) o„ a, md fa = ( pfcrO-i) on n, 



else [ else 

This coincides with the velocity basis in Example 1. If /c is varied in the coarse block, then 
ipij and ipik may not be orthogonal to each other any more. This can be observed from Figure 
m\ (right). 

3.2 Analysis of the mixed MsFV 

In this subsection we will show that the new mixed MsFV method is well defined. We have 
to check the conditions (J2]). The first one is straightforward since Vh = Qh and Uu = Vh- 
The second condition is verified below. We also give more details how to organize the 
computations. 

The mixed finite volume formulation ( 13. 8 p implies the following algebraic linear system 

AU + BP = , 

3.14 

CU — F, K ' 

where A is a mass matrix. Here B has only nonzero entries 1 and —1 and the BP represents 
the jump of P. The matrix C has only nonzero entries 1 and —1, and the sign depends 
upon the normal direction to which the corresponding flux entries of U associate. Because 
the first equation in ( 13. 8 p can be computed dual coarse block by dual coarse block and the 
support of velocity basis function lies in a dual coarse block, A can be represented as a block 
diagonal matrix, i.e., A = diag(A\, A 2 , ■ ■ ■ ), where each diagonal block Aj is the mass matrix 
associated to a dual coarse block. 

A straightforward calculation implies the following proposition. 



Proposition 3.2. Let B and C be defined in \3.1$ - Then B = C. 



Because the mass matrix A is block diagonal in the mixed MsFV, this allows to invert 
each block entry A\, A 2 , • • • to eliminate the flux U and the computation (for A^ 1 , A^ 1 , • • • ) 
is fast. It is known that mixed MsFEM also yields a system such as ( I3.14p . but mass matrix 
in mixed MsFEM is not block diagonal and one has to globally compute the inverse of the 
mass matrix to eliminate the flux U. In general, the computation for the inverse in mixed 
FEM is quite time-consuming when the mass matrix is large. This is an advantage of mixed 
MsFV over mixed MsFEM. Figure |3~41 shows the mass matrix for PTo mixed FEM (left) and 
mixed MsFEM (right), respectively, where 6 x 10 grid is used for both and permeability is 
heterogeneous (portion of SPE 10 in layer 85). 

We analyze the block diagonal entries of A. Let be the control volume with vertexes 
i, j, k, I. In DijM, u H can be represented by 

u h\d ijM = Uijipij + U ik tp ik + Ujitpji + U kl 1p k i. 
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Figure 3.4: Sparsity patterns for mass matrix in lowest order Raviart-Thomas mixed FEM 
(left) and sparsity patterns for mass matrix in mixed MsFEM (right) 

Consequently, the first equation in (13. 8p can be reduced to be in 
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(3.15) 



where a^^j := J D k ■ ipikdx and the other entries are defined similarly. We define A 
to be the most left matrix in (13. 15j) and it is a symmetric and positive Gram matrix. We 
would like to note that A is a representative for the diagonal entries in A defined in (I3.14p . 
We compute each multiscale basis functions ipij in fine scale by standard mixed FEM (e.g., 
Raviart-Thomas mixed FEM). Then ^ in D^u can be represented as iftij = J2 m rij t mipm, 
where ip^ are standard mixed finite element basis functions in fine scale. Hence A can be 
computed in the following way. 

Proposition 3.3. Let A h be the Gram matrix with entries (A h ) mn = f' /c -1 ^ • ip^dx, 
and each column of R consist of all rij )Tn . Then 

A = R T A h R. 

In particular, if permeability k is homogeneous (constant) in a dual coarse block and 
the mesh is uniformly square and no global information is used, then a straightforward 
calculation implies that the corresponding matrix 

A = diag(2, 2,2, 2), 

i.e., A is a diagonal matrix with diagonal entry 2. 

Because A is positive (or invertible), equation (13 . 1 5 j) implies that 



Uij 




~Pi 


-Pi 


u ik 


= A- X 


Pi 


-Pk 


Ujl 




Pi 


-Pi 


Ukl_ 




Pk 


-Pi_ 



(3.16) 



Here A 1 is a transmissibility matrix. We plug the expression (13. 16[) into the second equation 
in (13.81) to obtain a system about the pressure, 
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DP = F. (3.17) 

Here D = CA _1 B, where B and C are defined in (13.141) . We can show that D is sparse. In 
fact, from Equation (13. 8p . we know that only the fluxes of the interfaces of coarse volume 
contribute the mass conservation at volume Ki, and that each interface flux is determined 
by its neighbor pressures. Consequently, there are at most 9 nonzero entries for each row of 
matrix D for 2D rectangle cells. If we do not consider the restrictions of boundary condition, 
D is symmetric. This can be shown by using the fact B T = C. Our numerical studies show 
that D is positive. When permeability is homogeneous, a rigorous mathematics proof for 
the positiveness can be found in [25] . Hence equation (I3.17P is solvable. Once we obtain the 
pressure values, then we go back to equation (I3.16P to get the flux vector U. 

In particular, if the permeability k is a constant, the scheme in ( 13. 8 p coincides with the 
standard cell-centered difference scheme. In this case, we get the following matrix D for a 
2D uniformly square 3x3 grid, 
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Here we assume Neumann boundary condition without restriction j n pdx = 0. If we restrict 
J Q pdx = for a unique solution, e.g., replace the fifth row in D by (1, 1, 1, 1, 1, 1, 1, 1, 1), the 
matrix D becomes an invertible matrix D. 

We summarize the computation as following. 
Algorithm 1 

• For each interface Eij, we solve basis equation (13. lip . 

• By equation (13.81) . we formulate an algebraic system (I3.14p . 

• Eliminate U by local systems (13. 15)) and obtain equation (I3.17p . 

• Solve equation (13.171) to get pressure P and return to equation (13. 16)) to get U. 

• By basis equation (13.111) and U, downscale coarse scale velocity to fine scale velocity. 

3.3 Reconstruction of the fine-scale velocity field 

Fluxes across the interfaces of primary coarse grid can be accurately computed by (I3.16p . 
In many situations it is often needed to accurately compute the small-scale velocities (or 
fluxes) in regions of interest. Although we can obtain the velocity in fine grid by straight- 
forwardly prolonging the multiscale basis function into fine grid, the velocity are in general 
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discontinuous on the interfaces of the dual coarse grid. Then large errors can occur in the 
divergence field and local mass balance is violated. We can use a post-procedure to recon- 
struct a conservative fine-scale velocity which is continuous on the interfaces of both fine grid 
and coarse grid and fully consistent with the fluxes across the coarse grid interfaces given by 
the velocity multiscale basis functions. We would like to note that a similar post-procedure 
is used in MsFV [20]. We describe the post-procedure as following. 
Algorithm 2 

• Extract the velocities through the interfaces of all coarse primary cells K e Th- Let 
uji\dK be the velocity across dK, the boundary of K. 



• On each coarse grid K, the following local problem is solved by mixed finite element 
methods (or fine volume methods), 

f -div(fcVp) = / in K 

—kVp-n = UH\dK- n on 9K. 

• Define u H = @^ u^, where = —kVp and p solves (]3.19p . 

By the post-procedure, uh G H(div,Q), that is to say, Uh is continuous along all fine 
interfaces and coarse interfaces. 

Remark 3.3. Because there may exist some errors while computing uh, Jqk u h -ndx may be 
not close to j K fdx sufficiently We can replace the source term / in ( 13.191) by r^r J QK uu-ndx 
to remedy the little disparity for better accuracy. 

Remark 3.4. It the solution of problem (I3.19P is considered to be computationally expensive, 
another set of basis functions could be constructed only once and used during the simulation 
(see [20] for a similar approach and [27J for details). 

4 More mixed finite volume methods and their multi- 
scale analogues 

If we define different grids, approximation spaces or operators and div/j, then we can 
obtain different mixed finite volume methods. One mixed finite volume method has already 
been described in Example 1 of Section [21 and its multiscale analogue has been developed 
and analyzed in Section [3J In this section, we briefly present more examples of mixed finite 
volume methods and their multiscale analogues. 

The following Example 2 is closely related to Example 1. 

Example 2 (Mixed FV 2). The primary grid Th and the dual mesh T>h ore identical to the 
ones in Example 1. The discrete space Vh for pressure is the space of bilinear functions. The 
basis function of Vh for each cell center is one in the particular cell center and zero in all 
neighbors. Qh is the space of piecewise constants on the primary grid. The approximation 
space for the vector variable, Uh = Vh is the same as in Example 1 and Vh = V. 
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Figure 4.5: Voronoi box/Delaunay triangles 




Figure 4.6: Triangle D 

Unstructured grid is often used in practical simulations and mixed finite volume method 
can apply to the unstructured grid (see (13] for extensive discussions) . The following example 
describe a mixed finite volume method on an unstructured grid. 

Example 3 (Mixed FV 3) . We assume that the dual grid T>h consists of triangles and the cells 
in the primary mesh Th are control volumes around each vertex in T>h- Particular example 
is the Delauney mesh as a dual grid and the corresponding Voronoi grid as primary grid ( see 
Figure The space Vh is the space of piece- wise linear functions on the dual grid. The 

space Qh is the space of piece- wise constants on the primary grid. The approximation space 
for the vector variable, Uh = Vh, is the space of piece- wise constant vectors with continuous 
normal components and it is defined on the dual mesh (see Figure The construction 

of the basis is analogous with the procedure described in Example 1. The operator V/ t = V. 
Note that the dual grid can be any mesh of triangles, and the primary grid can consists of 
control volumes, not necessarily convex polygons, around each vertex. Frequently the control 
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Figure 4.7: Qaudrilateral mesh with volumes and covolumes 




Figure 4.8: Quadrilateral V{ 

volumes are formed by connecting the middle of the edges of each triangle with the center of 
mass of the triangle. The details of the method can be found in l2b^ . 

If the dual grid is Delauney mesh and the primary grid is Voronoi mesh, then we can 
approximate the scalar variable with piecewise constants, i.e., Vh is the space of piecewise 
constants on the primary grid. Then the operator Vh is defined by ( 12.61) . 

Here we consider a mixed finite volume method closely related to the multi-point flux 
approximation (MPFA) discretizations |24j . 

Example 4 (Mixed FV 4). The primary mesh is a general quadrilateral mesh (see Fig. 
\4- 7| j. The dual cells are formed by connecting the cell-centers with the centers of the edges 
of the quadrilaterals. These lines split each primary cell V\ into four quadrilaterals Vij, j = 
1, 2, 3, 4. The space Vh is the space of piecewise linear functions on with a common 
point in the cell center and the other two points on the edges on the quadrilateral V%j- (See 
Fig. \4-8j) . The space Qh is the space of piecewise constant functions on the primary grid. 
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The space Uh is the space of piecewise vector constants with continuous normals on the dual 
cells, and the space Vh is the space of piecewise vector constants on the dual cells. Note 
that these two spaces have different dimensions. Each local subspace Uh,i has four degrees of 
freedom. Each local subspace Vh,i has eight degrees of freedom. The operator Vh is defined 
by 

j ^hPh ■ v h dx = ^2 ^2 Vph ' Vh dx - 

Note that Vh here is not subspace of The approximation ofVp is nonconforming. 

One way to make the mixed finite volume method more robust is to use a conforming approx- 
imation. For example, it is possible to use as Vh the space of piecewise bilinear functions on 
Vij. One extra basis function is added for each intersection of four adjoint primary cells in 
two dimensions /2J/. Extra basis functions are required for three dimensional grids J77| /. 
We note that the same procedure works for the grids in Example 3 I3E/ . 

We used the mixed finite volume methodology to develop a new mixed multiscale method 
based on Example 1 in Section |3j We will sketch how we can follow the same procedure and 
derive mixed multiscale finite volume methods based on Example 2 - Example 4. Some of 
the methods derived following this framework are known, most are new. 

We assume that a consistent mixed finite volume method is defined on the fine mesh, i.e., 
fine primary and dual grids are given with the corresponding spaces and the operators. We 
need coarse primal and dual grids and we will assume that each coarse cell consists of fine cell 
of the same type, i.e., every coarse primary grid cell is a union of adjacent fine primary grid 
cells. The same is true for each coarse dual grid cell. This construction is straightforward for 
the structured grids. We suppose that an appropriate coarsening algorithm is used to define 
the coarse grid for the unstructured grids. The next step is to define the approximation 
spaces. We will require that the coarse discrete spaces provide some approximation of the 
corresponding functions. This requirement is easily fulfilled if we follow the same procedure 
on the coarse structured grids. For example the space Vh can consists of constant functions 
on the primary coarse grid Th- Then the space Uh has to be a multiscale space. We can 
consider a multiscale finite element (e.g., [HJ [181 ED] for detailed description) space for Vh 
and the space Uh could be the space of piecewise constant vector functions on the dual cells. 
The situation is more complicated for unstructured grids. We need approximation of the 
gradient of the pressure on unstructured grids and therefore we have to use a multiscale 
finite element space Vh- The space Uh can be either a standard piecewise vector constant 
space or a multiscale space. 

We provide below more details for several mixed multiscale methods that can be derived 
following our methodology and using the mixed finite volume methods on the fine grid in 
Examples 2 - Example 4. 

For better presentation, we call the mixed MsFV proposed in Section [3] to be Mixed 
MsFV 1. 

Example 5 (Mixed MsFV 2). The derivation below is related to Example 2 (Mixed MsFV 
2). The primary grid Th and the dual mesh Vh are identical to the ones in Mixed MsFV 
1 (Section^. We chose Vh to be the space of multiscale functions. The space Uh = Vh is 
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the space of piecewise constant vector functions with continuous normals on the dual coarse 
grid. The operator Vh = V. 

The basis ofVu on Th can be calculated in the same way as in f2^ . Then we will exactly 
reproduce the method proposed in JNftj . We can derive a different method by using a different 
basis for Vh- Tor example, the discrete harmonic basis J33\ [Lflj can be constructed in a 
multilevel way, that is cheaper, and the computations can be performed using only the fine 
matrix. 

Another method is derived if we select Vh to be the space of bilinear functions andUh = Vh 
to be a multiscale space defined in the same way as in Mixed MsFVl. The operator V h = V . 

The case for unstructured grids is more complicated. Here we sketch the derivation of a 
mixed MsFV for unstructured fine triangular grids. 

Example 6 (Mixed MsFV 3). We assume that the primary and dual grids Th, T>h are 
constructed from fine cells described in Example 3 (MFV 3) using an appropriate coarsening 
algorithm 133^ . The first method we propose uses a multiscale space Vh of discrete harmonic 
functions discussed in previous example. The space Uh = Vh is the space of constant vector 
functions with continuous normals. The operator Vh = V . Note that this method does not 
use global information. 

We recommend using the multiscale spaceUn = Vh when it is beneficial to use some global 
information about the problem and this information can be transferred using the velocity. 
The basis functions for Uh are constructed in a similar way as on the rectangular grid. We 
modify the definition of the velocity basis function A3.11\) as follows. Consider the triangle D 
on Figure \4~T\ The basis function corresponding to li is the solution of the following boundary 
value problem: 



' -div{kV(t)ij) = in A U Dj, 



kVifiij ■ n = < 



v(x)-n 
f v(x)-ndx 
v(x)-n 
f e v(x)-ndx 

jk 



(4.20) 








on l k U lj 
on e.ij U eji. 



Again Vh = V. 



Example 7 (Mixed MsFV 4). There are several ways to construct multiscale MPFA methods. 
One approach is presented in Wty . This requires a construction of a multisclale space Vh and 
the standard spaces forlAn and Vh- Another method with the ability to utilize the available 
global information can be constructed by using a multiscale space Vh, not necessarily the 
same as in IWf . and a multiscale spaces Uh and Vh defined in similar way as (I4.20p . There 
exists a third way for rectangular or quadrilateral grids, such that the edges of the coarse cells 
are straight lines. We can use the standard finite element space Vh and a multiscale space 
Uh- The framework also can be applied to the MPFA discretizations proposed in [77| 
and the corresponding mixed MsFV methods derived. 

We note that it is difficult to apply standard multiscale basis functions in production code 
because of the geometric grid information necessary to impose the appropriate boundary 
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conditions. It is also difficult to apply them for unstructured grids. If the coarse grid is 
very distorted, it may happen that a point where the four quadrilaterals meet is not in the 
support of the neighboring functions. This will decrease the accuracy. If discrete harmonic 
and multilevel basis is used, then computation becomes cheaper and can be applied on the 
matrix level (see [E]). The extensive study for these issues is still under investigation. 

5 Numerical results 

In this section, we apply the mixed MsFV propped in Section [3] to simulate incompressible 
two-phase flows in porous media. We will consider the porous media with non-separable 
scales and separable scales. In the first example, the permeability field is from SPE Compar- 
ative Solution Project [10] (also known as SPE 10) and its scales are non-separable, and it has 
channelized structure. In the second example, we consider the flows in two-point correlation 
permeability. The permeability is described with a two-point correlation function and has 
non-separable scales and distinct spatial variation. In the third example, the permeability is 
described by a periodic function with a small period and its scales are apparently separable. 
We apply the local mixed MsFV and global mixed MsFV to the flows in the three types of 
permeability fields. We will find that the mixed MsFV can provide accurate approximation 
on coarse grid and that using global information is able to greatly improve accuracy for the 
cases of non-separable scales. 

In our numerical simulations, we will perform two-phase flow and transport simulations. 
The equations are given (in the absence of gravity and capillary effects) by flow equations 

div(A(S)fcVp) = /, (5.21) 

where the total mobility \(S) is given by \(S) = \ W (S) + X (S) and / is a source term. 
Here, X W (S) = k rw (S)//i w and X (S) = k ro (S)/fi where \i D and \i w are viscosities of oil and 
water phases, correspondingly, and k rw (S) and k ro (S) are relative permeability of oil and 
water phases, correspondingly. The saturation is governed by 

— + div(F) = 0, (5.22) 

where F = vf w (S), with f w (S), the fractional flow of water, given by f w = X W /(X W + A G ), 
and the total velocity v by: 

v = Vw + v = -X(S)kVp. (5.23) 

In our simulations, we take k rw (S) = S 2 and k ro (S) = (1 — S) 2 . In the presence of capil- 
lary effects, an additional diffusion term is present in (15. 22ft and an efficient treatment of 
capillarity is proposed in [23]. We note that the porosity is 1 in the saturation equation 

We solve the two-phase flow system (15.211) and (15.221) by the classical IMPES (implicit 
pressure and explicit saturation). The saturation equation (I5.22p is discretized in fine grid 
by upwind finite volume method. The temporal discretisation is an implicit scheme, which 
is unconditionally stable but produce a nonlinear system (Newton-Raphson iteration solves 
the nonlinear system). For completeness, we describe the upwind finite volume method for 
equation (15.221) in Appendix [Al 
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We compare the saturation fields and water-cut data as a function of pore volume injected 
(PVI). The water-cut is defined as the fraction of water in the produced fluid and is given 
by q w /qt, where q t = q + q w , with q Q and q w being the flow rates of oil and water at the 
production edge of the model. In particular, q w — J anout f(S)v -nds, q t = J mout v -nds, where 
dVt out is the outer flow boundary. Pore volume injected, defined as PVI = y- f Q q t (r)dr, 
with V p being the total pore volume of the system, provides the dimensionless time for 
the displacement. We consider a traditional quarter five-spot problem, where the water is 
injected at left bottom corner and oil is produced at the right top corner of the rectangular 
domain. In all numerical simulations, multiscale basis functions are constructed once at the 
beginning of the computations. In the discussions, we refer to the grid where multiscale basis 
functions are constructed as a coarse grid. We use the global single-phase information (where 
X(S) = 1) to construct mixed MsFV basis functions. The global information is computed 
on fine grid at time zero and will not change throughout the simulation. 

In the simulations, we solve the pressure equation on the coarse grid by Algorithm 1 and 
use the post-procedure described in Algorithm 2 to re-construct the fine-scale velocity field 
which is used to solve the saturation equation. 

We solve the two-phase pressure equation f)5.2ip by the mixed MsFV (mixed FEM for 
reference solution). For the numerical simulations, we use 10 time steps for pressure equation, 
and for each pressure time step, we use 10 time steps to solve saturation equation. Hence, 
the time step for pressure is 0.1 PVI and the time step for saturation is 0.01 PVI. 

To assess the performance of the saturations and water-cuts obtained using the mixed 
MsFV, we compute the time-dependent pressure equation on fine grid by using lowest order 
Raviart-Thomas mixed finite element method, and this produces a reference velocity to solve 
a reference saturation solution S re f. By the reference saturation and the reference velocity, 
we get the reference water-cut W re f. We measure the relative saturation error in L 1 -norm 
and the relative water-cut error in L 2 -norm, 

\\SmsFV — SVe/IIWIISVe/IU 1 , H^AfsFy ~ W re f || L 2 / 1| W re f || £2 . 

where SmsFV and WmsFV denote the saturation and water-cut by the mixed MsFV, respec- 
tively. 

5.1 Flow in SPE 10 permeability 

For the first numerical example, we choose the SPE 10 permeability (layer 85), which is 
highly channelized and defined on a 60 x 220 find grid. The permeability map is depicted in 
Figure 15.91 We take 3x5 coarse grid for both the global mixed MsFV and the local mixed 
MsFV. We take tests for two different viscosity ratios of water and oil. 

For the first test of the example, we consider the case that viscosity ratio \i w j ' \i Q of water 
and oil is less than 1. Here we take fi w /fi = 1/10 for the simulation. Figure 15.101 shows the 
reference (fine-scale) saturation profile, the saturation profile using the global mixed MsFV 
and the saturation profile using the local mixed MsFV, respectively, at PV1=1. Figure 15.111 
shows the saturation error via different PVI times. From the figure, we find that the mixed 
MsFV using global information render more accurate saturation solutions throughout the 
whole simulation than the saturation from the local mixed MsFV. The water-cut curves are 



18 




Figure 5.9: Logarithm of SPE 10 permeability, layer 85. 

shown in Figure [5.121 for reference, global mixed MsFV and local mixed MsFV, respectively. 
The water break-through time is almost the same for the three methods, however, the water- 
cut curve by using global mixed MsFV is closer to the reference water-cut curve at early 
time than the local mixed MsFV. The average errors of saturation and water-cut are shown 
in Table [TJ We observe that the error of the solution of the global mixed MsFV is at least 
two times smaller than the error of the solution of the local MsFV. 

For the second test of the example, we consider the case when \i w j\i = 3. Figure I5TT31 
depicts the reference saturation profile and the saturation profiles uisng both the global mixed 
MsFV and the local mixed MsFV at PV1=1. Figure 15. 141 illustrates the saturation error via 
different PVI times and the water-cut curves as well. From Figure 15.131 and Figure 15.141 we 
find that global mixed MsFV performs better than local mixed MsFV in the point of view 
of accuracy and water-breakthrough time. Table [1] shows the average errors of saturation 
and water-cut for the test. 



Table 1: relative errors for saturation and water-cut, h w //j,q = 1/10 



mixed MsFV 


Water-Cut Error 


Saturation Error 


local mixed MsFV 


0.1140 


0.2024 


global mixed MsFV 


0.0510 


0.0870 



Table 2: relative errors for saturation and water-cut, fi w /fio = 3 



mixed MsFV 


Water-Cut Error 


Saturation Error 


local mixed MsFV 


0.1748 


0.2892 


global mixed MsFV 


0.1062 


0.1535 



5.2 Flow in two-point correlation permeability 

In the second example, we choose a realization of the permeability field generated using 
a two-point correlation function with correlation lengths in Xi-direction L\ = 0.4 and in 
^-direction L 2 = 0.05. Exponential variogram is selected (see e.g., [12]) to generate the 
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Figure 5.10: Saturation profiles, fi w /fJ>o = 1/10. Top: Reference saturation at PVF 
Middle: Saturation at PVI=1 by the global mixed MsFV. Bottom: Saturation at PVI=1 
the local mixed MsFV. 



relative saturation error 



- global mixed MsFV 



Figure 5.11: saturation error via time, fi w /fio = 1/10. 
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water-cut curve 



global mixed MsFV 
local mixed MsFV 




0.1 0.3 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

PVI lime 



Figure 5.12: water-cut curves, jiw/jio = 1/10. 



Reference saturation at PVI=1 




20 40 60 80 100 120 HO 160 180 200 



Saturation using global mixed MsFV at PVI=1 




20 40 60 80 100 120 140 160 180 200 220 



Saturation using local mixed MsFV at PVI=1 




20 40 60 80 100 120 140 160 180 200 220 



Figure 5.13: Saturation profiles, /i w / )i = 3. Top: Reference saturation at PVI=1 . Middle: 
Saturation at PVI=1 by the global mixed MsFV. Bottom: Saturation at PVI=1 by the local 
mixed MsFV. 
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Figure 5.14: water-cut curves and saturation error via time, fi w / fio = 3. 




Figure 5.15: Logarithm of two-point correlation permeability 

permeability. The permeability is defined on 200 x 200 fine-grid and depicted in Figure 15.151 
The viscosity ratio is n w /^ = 1/3 and the mixed MsFVs are implemented on 10 x 10 coarse 
grid. Figure 15.161 depicts the reference saturation, the saturation field using the global 
mixed MsFV and the saturation field using the local mixed MsFV, respectively. Figure 
15.171 demonstrate the relative saturation error via different PVI times. From the figure, we 
observe that mixed MsFV can provide a good approximation for the flow, and that using 
global information improves the accuracy. The water-cut curves are depicted in Figure 15.181 
for reference, global mixed MsFV and local mixed MsFV, respectively. Figure 15.181 shows 
that the water-cut curve in global mixed MsFV is more close to the reference water-cut curve 
than the mixed MsFV without using global information, and that the water break-through 
time in global mixed MsFV is almost the same as the reference water break-through time. 
The average errors of saturation and water-cut are shown in Table [3j From Table [31 we 
can observe: the saturation in the global mixed MsFV is almost 9 times better than the 
saturation in the local mixed MsFV, and the water-cut in global mixed MsFV is around 3 
times better than the water-cut in the local mixed MsFV. 
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Reference saturation at PVI=1 




Figure 5.16: Top: Reference saturation at PVI=1 for the second example. Middle: Satura- 
tion at PVI=1 by the global mixed MsFV. Bottom: Saturation at PVI=1 by the local mixed 
MsFV. 

relative saturation error 



Figure 5.17: saturation error via time for the second example 



Table 3: relative errors for saturation and water-cut for the second example 



mixed MsFV 


Water-Cut Error 


Saturation Error 


local mixed MsFV 


0.0464 


0.0610 


global mixed MsFV 


0.0127 


0.0075 
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Figure 5.18: water-cut curves for the second example 
5.3 Flow in periodic permeability 

In the third numerical example, we choose the permeability which is specified by a periodic 
function 

v ,yj 2 + 1.8sm(2vn//e) 2 + 1.8 cos(2vrx/e) ' 1 

The permeability is defined on 100 x 100 fine grid and its map is depicted in Figure 15.191 
Figure I5T2"U1 depicts the reference (fine-scale) saturation, the saturation using the global mixed 
MsFV and the saturation using the local mixed MsFV, respectively, at PVI=1. Here 5x5 
coarse grid is taken for simulation in both the global mixed MsFV and the local mixed MsFV 
and the viscosity ratio is ^ w /^i = 1/10. Figure 15.201 shows that the saturation profile by 
global mixed MsFV is almost the same as the saturation profile by local mixed MsFV and 
both of them are pretty close to the reference saturation profile. The saturation error via 
different PVI times is shown in Figure 15.211 It can be seen that the two saturation errors 
are very close and quite small. The water-cut curves are shown in Figure 15.221 for reference, 
global mixed MsFV and local mixed MsFV, respectively. Figure 15.221 shows that the three 
water-cut curves coincides each other. The average errors of saturation and water-cut are 
shown in Table HI We find from Table H] that: (1) the average errors of saturation and 
water-cut by global mixed MsFV and local mixed MsFV are comparable. (2) the average 
errors of saturation are less than 1% and the average errors of water-cut are less than 2%. 
This example shows that local mixed MsFV can provide very good accuracy in the case of 
separable scales (e.g., periodic) and its performance is as good as the global mixed MsFV. 
The example confirms the findings for many multiscale finite elements (e.g., [21 El l2Tj). 



Table 4: relative errors for saturation and water-cut for the third example 



mixed MsFV 


Water-Cut Error 


Saturation Error 


local mixed MsFV 


0.0074 


0.0171 


global mixed MsFV 


0.0082 


0.0169 
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Figure 5.19: Periodic permeability 



Reference saiuration at PVI=1 




10 20 30 40 50 60 70 SO 90 



Saturation using global mixed MsFV at PVI=1 




Saturation using local mixed MsFV at PVI=1 




10 20 30 40 50 60 70 SO 90 



Figure 5.20: Top: Reference saturation at PVI=1 for the third example. Middle: Saturation 
at PVI=1 by the global mixed MsFV. Bottom: Saturation at PVI=1 by the local mixed 
MsFV. 
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relative saturation error 



Figure 5.21: saturation error via time for the third example 



water-cut curve 
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0.2 - 
0.1 - 



Figure 5.22: water-cut curves for the third example 
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6 Conclusions 



We developed a framework of mixed MsFV methods for elliptic equations arising from flows 
in porous media. These methods take advantages of both mixed MsFEM and finite volume 
methods. The essential characteristic of mixed MsFv is the explicit approximation of both 
pressure and velocity. We proposed a new multiscale basis functions for velocity. Global 
information can be used to construct the multiscale velocity basis functions to improve ac- 
curacy for highly heterogeneous porous media. We analyze one of the mixed MsFV methods 
and apply it to simulate two-phase flows in porous media. We test the method on permeabil- 
ity fields with both non-separable and separable scales. Numerical examples demonstrate 
that the mixed MsFV can efficiently approximate the two-phase flows on coarse grid. Using 
global information in the mixed MsFV yields much better accuracy than the local mixed 
MsFV if the permeability field has strong no-local features. 

We also can use multiscale basis for the pressure. Discrete harmonic and energy mini- 
mizing basis constructed in a multilevel fashion is accurate and the computations are very 
efficient [33j E] • Moreover, the methods easily can be extended to unstructured grids with- 
out requiring extra geometric information. Further investigation of these issues is worth 
pursuing in the future. 



A Upwind finite volume method for Equation 15.22 



In the Appendix, we would like to present a finite volume discretization of the saturation 
equation (15.221) . Let 7^ be the common face (or edge) of iff (underlying fine grid) and Kj 
(underlying find grid) and n^- be the normal vector pointing from K% to Kj. Using the 
#-rule for temporal discretization and a finite- volume scheme for the saturation equation, it 
follows the following form. 



i_(S r +i _ s») + _L J2[9 FlJ (S n+1 ) + (1 - e)F ii {S n )] = 0, (A.24) 
where S™ is the cell-average of water saturation at t — t n , i.e., 

S™ = {S[X,t n )) K h 

and Fij is a numerical approximation of the flux over 7^, i.e., 



Fij(S) « / f w {S)ijU i: j ■ riijds. 



Here f w (S)ij denotes the fractional-flow function associated with 7^ and the first-order 
upstream weighting scheme for it is defined as 

f /m / fw(Si) if U ■ Tlij > 

M )ij \ f w (Sj) if 0. 
For 6 = or 1, we can write (1A.24P vector form 

S n+1 = S n + (S t x fAf(S m ), m = n or n + 1, 
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where = 

If 9 = 0, then ( 1A.24}) is an explicit scheme and only stable provided that time step At 
satisfies a stability condition. 

For 9 = 1, (1A.24|) is an implicit scheme and unconditionally stable but gives rise to a 
nonlinear system. Such a nonlinear system is often solved with a Newton-Raphson iterative 
method. Define 

G(S n+1 ) = S n+1 -S n - (S f x ) T Af(S n+1 ). (A.25) 
By Taylor expansion, we have 

G(S n+1 ) « G(S n ) + G'{S n ){S n+1 - S n ). 

Noticing G(S n+1 ) = we have 5S n := S n+1 - S n = -[G'(S n )}^G(S n ). From (TA~25|) . we 
have 

C(5) = J-(^) T Af # (5), 

where f'(S)i = f'(Si). Hence 

S n+i = S n + 5S n . 

This iteration proceeds until pre-defined iterations are reached or the norm of 5S n is smaller 
than a prescribed value. 
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